Q5: Planetary orbits

We want to consider planetary orbits. To do this, we need to solve Newton’s second law together with Newton’s law of gravity. If we restrict ourselves to the x-y plane, then there are 4 quantities we need to solve for: \(x\), \(y\), \(v_x\), and \(v_y\). These evolve according to:

\begin{align*} \frac{dx}{dt} &= v_x \\ \frac{dy}{dt} &= v_y \\ \frac{dv_x}{dt} &= a_x = -\frac{GM_\star x}{r^3} \\ \frac{dv_y}{dt} &= a_y = -\frac{GM_\star y}{r^3} \end{align*}

To integrate these forward in time, we need an initial condition for each quantity. We’ll setup our system such that the Sun is at the origin (that will be one focus), and the planet begins at perihelion and orbits counterclockwise.

The distance of perihelion from the focus is:

\[r_p = a (1 - e)\]

where \(a\) is the semi-major axis and \(e\) is the eccentricity. The perihelion velocity is all in the \(y\) direction and is:

\[v_y = v_p = \sqrt{\frac{GM_\star}{a} \frac{1+e}{1-e}}\]

We’ll work in units of AU, years, and solar masses, in which case, \(GM_\star = 4\pi^2\) (for the Sun).

Your initial conditions should be:

• \(x(t=0) = r_p\)

• \(y(t=0) = 0\)

• \(v_x(t=0) = 0\)

• \(v_y(t=0) = v_p\)

Use the scipy ODE integration methods to integrate an orbit and plot it

from BicoccaCoursePython2024.quarta_lezione import planet

Earth = planet()

Earth.integrate_motion_equations()

ani= Earth.animate_jnb()

ani

Once Loop Reflect

Q7: Noisy signal

A convolution is defined as:

\[(f \star g)(t) \equiv \int_{-\infty}^{\infty} f(\tau) g(t - \tau) d\tau \tag{1}\]

It is easy to compute this with FFTs, via the convolution theorem,

\[\mathcal{F}\{f \star g\} = \mathcal{F}\{f\} \, \mathcal{F}\{g\} \tag{2}\]

That is: the Fourier transform of the convolution of \(f\) and \(g\) is simply the product of the individual transforms of \(f\) and \(g\). This allows us to compute the convolution via multiplication in Fourier space and then take the inverse transform, \(\mathcal{F}^{-1}\{\}\), to recover the convolution in real space:

\[f \star g = \mathcal{F}^{-1}\{ \mathcal{F}\{f\} \, \mathcal{F}\{g\}\} \tag{3}\]

A common use of a convolution is to smooth noisy data, for example by convolving noisy data with a Gaussian. We’ll do that here.

Here’s some noisy data we’ll work with

import numpy as np
import matplotlib.pyplot as plt

def fdata(x, L):
    A = L/10.0
    return 2*np.sin(2*np.pi*x/L) + x*(L-x)**2/L**3 * np.cos(x) + \
           5*x*(L-x)/L**2 + A/2 + 0.1*A*np.sin(13*np.pi*x/L)

N = 2048
L = 50.0
x = np.linspace(0, L, N, endpoint=False)
orig = fdata(x, L)
noisy = orig + 0.5*np.random.randn(N)

plt.plot(x, noisy)
plt.plot(x, orig)

[<matplotlib.lines.Line2D at 0x2204da69790>]

SciPy provides a convolution function scipy.signal.convolve() that can do the convolution for us directly. To smooth the data, we want to use a Gaussian, which can be produced by scipy.signal.gaussian().

Convolve the noisy data with a Gaussian and plot the result together with the original data orig. You’ll need to play with the width of the Gaussian to get a nice smoothing. You also will need to normalize the Gaussian so that it sums to 1, otherwise, your convolved data will be shifted verfically from the original function.

from BicoccaCoursePython2024.quarta_lezione import signal
sig= signal(noisy,x)
sig.clean()
_,ax=sig.plot()
ax.plot(x,orig,label='Original function',linestyle='-.')

[<matplotlib.lines.Line2D at 0x220497cd650>]